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Systems with long-range interactions show a variety of intriguing properties: they typically ac- 
commodate many meta-stable states, they can give rise to spontaneous formation of supersolids, 
and they can lead to counterintuitive thermodynamic behavior. However, the increased complexity 
that comes with long-range interactions strongly hinders theoretical studies. This makes a quantum 
simulator for long-range models highly desirable. Here, we show that a chain of trapped ions can be 
used to quantum simulate a one-dimensional model of hard-core bosons with dipolar off-site interac- 
tion and tunneling, equivalent to a dipolar XXZ spin- 1/2 chain. We explore the rich phase diagram 
of this model in detail, employing perturbative mean-field theory, exact diagonalization, and quasi- 
exact numerical techniques (density-matrix renormalization group and infinite time evolving block 
decimation). We find that the complete devil's staircase — an infinite sequence of crystal states 
existing at vanishing tunneling — spreads to a succession of lobes similar to the Mott-lobes found in 
Bose-Hubbard models. Investigating the melting of these crystal states at increased tunneling, we 
do not find (contrary to similar two-dimensional models) clear indications of supersolid behavior in 
the region around the melting transition. However, we find that inside the insulating lobes there are 
quasi-long range (algebraic) correlations, opposed to models with nearest-neighbor tunneling which 
show exponential decay of correlations. 

PACS numbers: 75.10.Jm, 03.67.Ac, 67.85.-d, 02.70.-c 



I. INTRODUCTION 

Dipolar interactions in a lattice system can lead to 
exciting new physics not present in systems with short- 
range interactions [Q. For example, they can give rise 
to new ground states, like crystal states with fractional 
filling factors or supersolids. For instance, on the square 
lattice under hole-doping from half- filling, dipolar inter- 
actions can stabilize a supersolid phase [2], in contrast 
to interactions of limited range [3H5]. Further, dipolar 
interactions can lead to the appearance of a large num- 
ber of metastable states [5J. On the one hand, these can 
complicate finding the ground state in numerical studies, 
while on the other hand they might be useful for appli- 
cations in quantum information and as quantum mem- 
ories [7H9]. Systems with strong long-range interactions 
(where the integral over the interactions does not con- 
verge) can even exhibit counterintuitive thermodynamic 
effects like super-extensive behavior, breaking of ergodic- 
ity, or a non-concave entropy [TU]. These effects, however, 
are not the subject of the present work, since we consider 
weak long-range interactions. 

Since long-range interactions can make theoretical cal- 
culations more difficult, it would be highly desirable to 
study all this extremely rich physics with an analog quan- 
tum simulator - specially for dynamics or two or more 
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dimensions. In one dimension (which is the subject of 
this work) , standard analytical and numerical techniques 
typically suffice to give a good picture of the phase di- 
agram, but this is an important limit to gain intuition. 
Its experimental implementation allows to judge the reli- 
ability and performance of analogue quantum simulators. 
Natural candidate systems for this type of quantum sim- 
ulation are trapped ions, where the high degree of control 
over state preparation, evolution, and readout allows ac- 
cess to a wide range of models. In 2004, trapped ions 
were first proposed for quantum simulations of certain 
classes of spin models by Porras and Cirac |TT] . In a dif- 
ferent context, similar spin models were derived earlier 
by Mintert and Wunderlich [T^], who studied the possi- 
bility of individual ion addressing using inhomogeneous 
magnetic fields. Recent proof-of-principle experiments 
showed that, indeed, quantum spin systems can be sim- 
ulated faithfully by systems of cold trapped ions [131 [Tl] . 
For a review of trapped ion quantum simulations see [TS] . 

An advantage of trapped ions is that dipolar interac- 
tions (which are the most common type of long-range 
interactions) are naturally achieved without additional 
experimental effort since the interactions are medi- 
ated by phonons. This leads also to the peculiarity that, 
contrary to other possible setups for the quantum simula- 
tion of spin models, not only two-body (density-density) 
interactions, but also tunneling terms can be made long 
ranged. 

In this work, we consider an experimentally feasible 
setup which allows to quantum simulate a chain of XXZ 
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spins with dipolar interactions in a magnetic field. The 
system, which, as we show, supports a rich ground-state 
phase diagram, is described by the Hamiltonian 

H = JJ2 -T^T3 [cos eSfSf + an6(S?S? + S^S/)] 

(i) 

where fi is the chemical potential (equivalent to an ex- 
ternal magnetic field), and the Sf are spin-operators at 
site i. 

The problem can be mapped to a system of hard-core 
bosons using the standard Holstein-Primakoff transfor- 
mation, Sf — > rii — 1/2, Sf — > a], and S~ — > a i; where 
a\ (cbi) creates (destroys) a boson at site i and where n, 
denotes the corresponding number operator. In this pic- 
ture, an up-spin corresponds to an occupied site and a 
down-spin to an empty one, the ZZ terms translate to 
dipolar off-site density-density interactions, and the XX 
and YY interactions become long-range tunneling terms. 
Since typically tunneling terms are short-ranged, the lat- 
ter has to be seen as an important peculiarity of our 
model. In the following, we will use both the spin and 
boson pictures interchangeably, because some aspects are 
better described in terms of spins, while others are more 
familiar in the language of hard-core bosons. 

By varying the angle 9, we can explore all ranges of 
relative strength of the interactions, e.g. the XX chain 
in a transverse field (9 — tt/2), or the Ising model (9 = 
0) — but both with long-range interaction. Moreover, 
the model can be tuned from negative to positive XY 
interaction, with the latter leading to a deformation of 
the phase diagram due to frustration effects. 

At 9 — 0, where Hamiltonian ([lj behaves essentially 
classical, this model is known to show a succession of in- 
sulating states with different rational filling factors |16j . 
This succession, called 'devil's staircase', covers the en- 
tire range of the chemical potential, with (for an infinite 
system) each rational filling factor occuring on a finite 
range of chemical potential. Since the Hamiltonian ([I]) 
shows a symmetry of up-down spins (or in hard-core bo- 
son terminology, a particle-hole symmetry), we consider 
in this work only negative magnetization (filling factors 
lower than 1/2). When \9\ grows, the insulating states be- 
come destabilized and finally melt to a superfluid phase. 
We concentrate on characterizing the complete phase di- 
agram of the model, and the nature of the correlations in 
each phase. We find that the long-range nature of the in- 
teractions changes qualitatively one important aspect of 
the physics of the problem: inside the insulating phases, 
the off-diagonal correlations decay algebraically instead 
of exponentially (as is typical in models with short-range 
tunneling) [17 . We call this phase a quasi-supersolid, 
since this is the closest analogue to a supersolid that can 
exist in one dimension: A supersolid is characterized by 
the coexistence of diagonal and off-diagonal long-range 
order (LRO). In ID, however, off-diagonal LRO cannot 



occur [18J, and the 'best' that can exist is the coexistence 
of diagonal LRO and off-diagonal quasi-LRO with alge- 
braically decaying correlations. This quasi-supersolid is 
exceptional, since normally systems with dipolar interac- 
tions in ID can be described by Luttinger theory |19l |2T)] . 
In that case, diagonal correlations decay with the Lut- 
tinger parameter K, and off-diagonal correlations with 
1/K. This means that if one of them decays slowly, the 
other one decays very fast. In contrast, in our model with 
long-range tunneling both diagonal and off-diagonal cor- 
relations decay slowly, and Luttinger theory is not a valid 
description. 

Recently, a similar model with short-range tunneling 
has been investigated by a strong-coupling expansion in 
the context of ultracold dipolar atomic gases |21| . Placed 
in a one-dimensional optical lattice, in the limit of strong 
on-site interaction, such a system becomes similar to 
the one described by Hamiltonian 0, but with nearest- 
neighbor (NN) instead of long-ranged tunneling. Hence, 
the quasi-supersolid phase cannot be observed in this sys- 
tem. Throughout the paper, we draw comparisons be- 
tween our model and the model with NN tunneling and 
dipolar interaction, as well as with a model with both 
NN tunneling and interaction (the NN-XXZ model). 

We organize this paper as follows: First, we briefly 
explain in Sec. [TT] how the desired model Hamiltonian H 
can be implemented in a system of trapped ions. Then, in 
Sec. |III| we outline our expectations by discussing previ- 
ous related results. The following sections are dedicated 
to a thorough investigation of the ground state phase di- 
agram of the model. Our most accurate analysis of the 
problem comes from numerical density matrix renormal- 
ization group (DMRG), Sec. [V] However, it is good to 
first form intuition via analytical approaches, even if ap- 
proximate. We will present mean field and pcrturbative 
results in Sec. |IV| to find the upper borders of the crys- 
tal lobes. Some limitations of DMRG can be overcome 
with other numerical techniques. For instance, we can 
study infinite systems (as opposed to finite sized) with 
the infinite time evolving block decimation (iTEBD) al- 
gorithm (section Sec. VI). Finally, experimentally rele- 
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vant small systems can be thoroughly studied with exact 
diagonalization (Sec. 
final remarks. 



II. EXPERIMENTAL IMPLEMENTATION 

In this section, we shortly explain how the desired 
Hamiltonian, Eq. 0, can be simulated in an ion trap 
experiment. The discussion is a slight generalization of 
Ref. [TT], to which we refer the reader for details. By ap- 
plying a standing laser wave to a chain of trapped ions, 
one can arrive (after a canonical transformation) at an 
effective spin model where two internal hyper-fine states 
act as a pseudo-spin. In this scenario, the Coulomb in- 
teraction between two ions i and j transmits an effective 
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spin-spin interaction 



j a — 



a«. ( j_ 

m \IC a 



(2) 



where F a is the field strength of the laser propagating 
along direction a = x,y,z, m is the ion mass, and the 
elasticity matrix is 



K H = 



Ee 2 / rn 



+c a - 



2 / 

e jm 



* = j 



(3) 



Here c^.y = 1, c z = —2, and oj a is the vibration fre- 
quency of an individual ion. In the limit where the 
energy scale given by the Coulomb interaction between 
neighboring ions becomes much smaller than the en- 
ergy scale of the vibration of the individual ions u) a , i.e. 
f3 a = \c a \ e 2 /mw^dj <C 1, the decay of the spin-spin in- 
teractions obeys a dipolar power law. Here, we intro- 
duced the mean inter-ion distance do, and assumed that 
the ion chain is homogeneous. In a realistic experimen- 
tal situation, if an overall trapping potential is used, the 
inter-ion spacing becomes non-uniform. In the central 
region of a long ion chain, however, the inhomogeneity 
can be neglected. Another way of creating an equidis- 
tant chain would consist of placing the ions in individual 
microtraps [22 . For simplicity, we will focus on the ho- 
mogeneous situation. 

Using three pairs of laser waves, one can thus induce 
the desired dipolar spin interactions in all the three di- 
rections of spin space. Choosing the detuning of the laser 
beams appropriately, one can tune the spin-spin coupling 
to negative or to positive, as desired. The missing chem- 
ical potential term can be achieved rather easily via an 
rf-frequency field driving the transition between the two 
hyper-fine states which consitute the effective spin 1/2. 
This completes all terms in Hamiltonian (JT|). 



III. EXPECTED BEHAVIOR OF THE MODEL 

Before proceeding to the detailed numerical analysis of 
the ground state phase diagram, let us briefly sketch the 
expected behavior of the present model, starting from 
previous related results. 

For 9 = 0, Hamiltonian Q describes a ID system of 
localized hard-core bosons with repulsive dipolar inter- 
action, a classical model that can be solved analytically 
|16j . Due to the long-range nature of the interactions, 
for any given filling factor, the particles arrange in a 
periodic crystal pattern (a generalized Wigner lattice). 
For a given filling factor, these periodic patterns can be 
constructed by maximizing the mutual distance between 
the particles. Every rational filling factor q = — occu- 
pies a finite extent in fi, thus giving rise to a plateau of 
fixed particle density, and these plateaux cover the entire 
range of /i. Plotting the filling factor against the chemical 



potential /i yields a self-similar structure — a complete 
devil's staircase, similar to Ref. [16_. The stair's steps 
have a width 
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This means that fillings with a large denominator - 
equivalent to a large crystal period — become the ground 
state only in a very small range of chemical potential. 

A finite tunneling, i.e. (0mod7r) ^ 0, introduces quan- 
tum mechanics into the system. The tunneling allows 
particles to gain kinetic energy, which destabilizes the 
crystal until it melts at some critical tunneling strength. 
This gives rise to a superfluid (SF) phase where the par- 
ticles are delocalized over the chain. 

There are two possible scenarios for a transition from 
the crystal structure to the molten phase: a direct 
crystal-SF transition, or an intervening supersolid phase. 
The latter is an exotic quantum phase showing crystal 
and SF behavior at the same time, which has been pre- 
dicted for similar — but two-dimensional — hard-core 
boson systems [5j |2"3"H2"T] . Until recently, there was only 
one claim of an experimental realization of a supersolid 
gum] in 4 He, which is still disputed (SHI EH- Therefore, 
it would be interesting to find other systems which show 
supersolid behavior and are cleaner to interpret. One 
such experiment has been carried out very recently |32j . 
where an atom cloud in a cavity breaks spatial symmetry 
while at the same time showing off-diagonal long-range 
order. 

In the limit of 9 = —tt/2, Hamiltonian ([I]) describes 
a ferromagnetic XY model. Here, the most important 
effects of the long-range nature of the tunneling can be 
captured in a renormalized nearest-neighbor (NN) inter- 
action: first, all interactions work towards the same or- 
dering, and, second, in one dimension the integral over 
dipolar interactions converges. For 8 = +tt/2 the system 
is a frustrated XY antiferromagnet and the behavior is 
less obvious. It turns out, however, that a similar renor- 
malization to a NN interaction captures the main physics. 
Hence, analogies to a NN-XY model can help understand- 
ing the behavior of the system near 9 = ±7r/2. 

In the following, we study how the ground state phase 
diagram of Hamiltonian ([TJ develops function of 9. 
We start with two approximative methods to obtain qual- 
itative insights. 



IV. MEAN-FIELD APPROXIMATIONS 

A. Perturbative mean-field theory 

A first rough understanding of the crystal-SF transi- 
tion can be obtained by a perturbative mean-field theory 
(PMFT), valid for small tunneling. While not being very 
accurate in ID, such a mean-field treatment generally 
becomes better for longer-ranged interactions. 
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Following the standard recipe of this Ansatz, we de- 
compose the Hamiltonian into two parts, H = Ho + Hi. 
We choose 



H = jJ^—^SfS?- M 
i<j \*~3 



3 1 3 cos (6) 



XX, ( 5 ) 



i.e. we consider the insulating states, the exact ground 
states at the devil's staircase (6 — 0), as the unperturbed 
states. Afterwards, we introduce the perturbation 
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H x = J tern 6 J2 X ., (S?S? + S?SP 



(G) 



= J tan 9 V — L_^5+Sfr. 

The aim is now to calculate the expectation value of the 
SF order parameter 



<Pi = (Si > = Tr (pS- 



(7) 



where S~ is the spin lowering operator at site i and p 
is the density matrix of the system. At the tunneling 
strength where tpi becomes finite, the assumption that 
the ground state is localized is no longer valid. This 
gives an upper border for the crystal lobes. 

For the evaluation of the expectation value in Eq. 
we need the density matrix p. which is given by 



P 



- ~PH 



/z, 



(8) 



with Z the partition function and (3 — 1/ (UbT) the 
inverse temperature. We are only interested in ground 
state properties, equivalent to the limit (3 — > oo. In this 
limit, Z — > e _/3E °, where Eq is the ground state energy. 
For small tunneling, i.e. small Hi, we can approximate 
p using the Dyson expansion 



-f3H 



-pH 



-PH 



P 



dTe THo Hie 



-tH 



(9) 



In mean-field approximation, we have 
HicJ2 (S+ (Si) + (Si) S- (Sr)) 



i»-jT 



(10) 

Inserting Eqs. ^ to ( 10 1 into we obtain the following 
self-consistent equations for the SF order parameter 



I Jtan0 v ^ 

x ^- L S t*i 



where we introduced the abbreviation 
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-tH 



,(n) 



(12) 



It can be checked that in the trace only one-hole excited 
states and the ground state itself have to be taken into 
account. For small J tan 9, Eq. (Ill has only the trivial 
solution tfi = 0. At some critical J c (p), however, SF 
order may develop, i.e. ifi =/= 0. This is equivalent to 
an instability under adding or removing a single particle. 
The point where this happens gives an upper border for 
the crystal phase. The main disadvantage of the current 
method is that there are potentially more complicated 
excitations, e.g. the addition of a particle plus a reloca- 
tion of the neighboring particles. Such a deformation of 
the crystal could decrease the potential energy. In the 
current simple Ansatz, however, this type of excitations 
is not captured. 

For our calculations, we assumed an infinite system 
with the restriction that the states have a periodicity of 
12 sites. Since the crystal phases with a large periodic- 
ity become very small in their extent in p [by virtue of 
Eq. Q], this is a reasonable restriction which still cap- 
tures the most prominent features of the phase diagram. 

The insulating-lobe structure that we obtain is shown 
in Fig. [T] The thick line follows the breakdown of the 
ground state. A dotted line marks the values of p where 
the ground state magnetization changes, i.e. where a dif- 
ferent crystal structure becomes lower in energy. How- 
ever, even when a given crystal structure ceases to be 
the ground state, it can still be a metastable state with 
respect to one-particle or one-hole excitations. In Fig.[T] 
thin lines mark the corresponding regions where the most 
important crystal states are metastable. For clarity, 
states that constitute the ground state over only a very 
small region of the phase diagram have been excluded. 
Figure [T] shows the rich structure of metastable crystal 
states. Recently — in the context of ultra-cold dipo- 
lar neutral atoms — it has been proposed to use such 
metastable states as quantum memories [H1IH]- Another 
feature of the crystal structure is that it is more stable for 
frustrated antiferromagnetic XY-interaction (9>0) than 
for ferromagnetic XY-interaction (9<0). 

We also computed the phase diagram in a Gutzwiller 
mean-field theory (not shown) . The qualitative behavior 
is similar, but the lobes are somewhat smaller. The rea- 
son is that PMFT considers destabilization under single- 
particle or single-hole excitations, while Gutzwiller mean- 
field theory captures better more complicated excita- 
tions. 



B. Wigner-crystal melting at low filling 

For very low magnetizations, we can draw an anal- 
ogy to the melting of Wigner crystals |3.3I [34] . As de- 
scribed in Sec. |III| at 9 = 0, the hard-core bosons are 
perfectly localized, with an inter-particle distance given 
by the filling fraction. At finite tunneling, 9 i 0, the par- 
ticles spread, but at small tunneling and low filling it is 
reasonable to assume that they remain spatially well sep- 
arated. Under this assumption, we can approximate the 
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Figure 1: Stability regions of the crystal states within PMFT. 
The thick line traces the melting of the crystal ground state. 
The dotted lines mark the values of [i where the crystal ground 
state changes periodicity. The thin lines delimit the regions 
where the crystal states are metastable. Note the asymmetry 
between the negative (ferromagnetic) and positive (antiferro- 
magnetic) 8 side. F-SF stands for ferromagnetic and AF-SF 
for antiferromagnetic superfiuid. 




Figure 2: The Lindemann parameter (described in the text) 
as a function of the interaction parameter 9 and the mean 
distance between particles in the crystalline phase. The white 
lines are contour lines at L — 0.02 ... 0.14 in steps of 0.02. For 
comparison two black lines following the law tan# oc d^ 1 are 
also shown. 



total wave-function by a product of Gaussians represent- 
ing the individual particles |43| . Using this Ansatz, one 
can obtain a self-consistent equation for the distribution 
of any particle over the lattice. 

The spread of the wave packets increases with 8 until 
the initial assumption that the particles are well sepa- 
rated breaks down, and the crystal has to be considered 
as molten. A me asure for thi s is the normalized variance 

L = ^ — jj^ J (x 2 ) — (x) 2 , called the Lindemann pa- 
rameter |35| . As shown in Fig. [2] for large inter-particle 
spacing d 0l the system behaves similar to the continuum 
limit, where it is known [33 that there is a melting tran- 
sition at Jtan# c = const /do. We find that the contour 
lines of the Lindemann parameter from our crude approx- 
imation follow this behavior well for large <i - At small 
inter-particle distances there are deviations, but here the 
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Figure 3: Mean polarization (Z) /N = \ J2f=i (°j ) l N for a 
system with dipolar hopping and TV = 60. The log-log plots 
in Panels A to D show, for N = 102, the decay of the in-plane 
correlations G (j) at the marked points of the phase diagram, 
comparing dipolar hopping (points) with NN hopping (solid 
lines). The dashed lines are guides to eye proportional to j -3 . 
For dipolar hopping, the decay is algebraic everywhere, even 
within the lobes (A and C), where it follows the dipolar ex- 
ponent a = 3, while for NN hopping the decay is exponential. 
In the SF (B and D), the algebraic decay is faster than for 
NN hopping on the antiferromagnetic side and slower on the 
ferromagnetic side. Panel E shows a zoom on the dashed sec- 
tion together with the finite-size scaling result of the crystal 
lobe at 1/3 filling, as explained in the text. 



assumptions made in our approximation are not valid 
anyway. The physical explanation behind the observed 
behavior is that for large inter-particle distances the re- 
pulsive interaction becomes very weak. Hence, for larger 
d a smaller gain in kinetic energy suffices to destabilize 
the crystal. 



V. DENSITY-MATRIX RENORMALIZATION 
GROUP (DMRG) 

While mean-field theories can provide some physical 
understanding of the system, in ID they are not very 
precise. Therefore, we turn now to a thorough numerical 
analysis of the phase diagram by the quasi-exact density- 
matrix renormalization group (DMRG) method [3SH3S]. 
We consider chains with open boundary conditions of up 
to N = 102 sites, and with a bond dimension of D = 
128. The interaction range was always cut-off at half the 
length of the chain. 

Figure p] presents the mean polarization (Z) /N = 
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cr|) /JV, where AT is the number of sites and a* 



the Pauli-z matrix of spin j. We find that for large 
enough field fj,, or strong ferromagnetic ZZ-interaction 
(corresponding to 9 not too far from ±7r), the system is 
in a fully polarized state, quantitatively similar to what 
we saw in the mean-field calculations. However, as ex- 
pected, the general precision of the mean-field calcula- 
tions in this one-dimensional system is low: the DMRG 
results indicate that the #-range of the crystal lobes is up 
to an order of magnitude smaller than predicted in mean- 
field. In fact, in the global view given in the main panel 
of Fig. [3j only the plateau for 1/2 filling is discernible. 
Panel E of Fig. [3] shows the small region of the phase di- 
agram where the 1/3-fiHing crystal lobe is located. Due 
to the long-range nature of the hopping terms, the po- 
larization is asymmetric with respect to ±9, in contrast 
to [ST]- 

We find that open boundary conditions play an impor- 
tant role for the finite systems used here. For example, 
in the center of the main panel of Fig. [3] there appears a 
broad plateau with (1/2— 1/N) filling (i.e. with one spin 
flipped away from half filling). In the thermodynamic 
limit, this phase is not as stable and merges with the 1/2- 
filling plateau to a single lobe. This peculiarity is highly 
relevant for the correct interpretation of future experi- 
ments, which might be carried out with open boundary 
conditions for technical convenience. Moreover, bound- 
ary effects play an important role in the determination of 
the decay of correlations, which we now turn to describe. 

The log-log plots in panels A to D show the decay of 

the in-plane correlations G (j) = ( f7 jv/ 2 cr JV/2+j) a ^ some 
selected points in the phase diagram (filled circles), and 
compare the result to the same system but with nearest- 
neighbor (NN) hopping (solid line). In the superfluid 
phase (panels B and D), the decay is algebraic for both 
NN and dipolar hopping. Panels A and C lie inside the 
lobes corresponding to 1/2 and 1/3 filling, respectively, 
which is visible through the oscillation of the correlation 
function superposed with an overall decay. We see that, 
for dipolar hopping, the decay is algebraic everywhere, 
even within the lobes, where for NN hopping the cor- 
relations decay exponentially. In the lobes, the decay 
for dipolar hopping follows the interactions with expo- 
nent a = 3 [17] ■ Notice that a clear fit for this expo- 
nent can only be obtained for rather large chains with 
A^ > 60 spins. Even more, although the strength of the 
XY-interactions decays very rapidly with distance, the 
power-law decay of correlations inside the lobes cannot 
be observed at all if the interaction is truncated. 

Since inside the lobes we have a crystal structure, 
the long-range decay of the transverse correlations is an 
anomaly when compared to other models with insulating 
phases. For this similarity to a supersolid, we call this 
phase a quasi-supersolid phase. In a true supersolid, a 
crystal structure coexists with a long-range order in the 
transverse correlations. However, in ID off-diagonal LRO 
cannot be spontaneously broken [13]. Hence, the closest 
analogue to a supersolid that can exist in ID is a state as 




Figure 4: Exponents a of an algebraic fit G a (j) = coj a to 
the correlations G(j), for A^ = 60. The panels show cuts at 
fj. = (1/2 filling, A) and 9 = tt/2 (XY model, B). The black 
lines in A are the exact results for the NN XXZ model. Anti- 
ferromagnetic interaction induces larger exponents compared 
to the NN XXZ model, while ferromagnetic interaction leads 
to smaller exponents. Inside the insulating lobes, a — 3 (devi- 
ations are due to finite size effects). Results for NN tunneling 
are shown in red (dashed line in A, filled circles in B). 



the present one — with diagonal LRO and algebraically 
decaying off-diagonal correlations. 

A change in the behavior of the correlations marks the 
transition from a crystal state to the SF. For a quanti- 
tative evaluation, at each point in the phase diagram we 
fit an algebraic trial function G a (j) = coj~ a . The value 
of a for a chain of N = 60 spins is shown in Fig. |4j The 
cuts of Fig. gat n = (1/2 filling, A) and 9 = tt/2 (XY 
model, B) show the wide range that the exponent of the 
correlations takes in this system. They also demonstrate 
the influence of frustration: the correlations decay faster 
at the antiferromagnetic side (9 > 0) than at the ferro- 
magnetic side (9 < 0). The exponents are similar to the 
NN XXZ model (black line in A), but cover a broader 
range. In particular, they are not independent of the fill- 
ing for 9 = n/2 (panel B), where the NN XXZ model 
has a = —0.5. The results for a system with dipolar 
ZZ-interaction but NN XY-interaction are shown in red 
(dashed line in A, filled circles in B). The main qualita- 
tive difference is that inside the insulating lobes the decay 
is no longer a power law, but follows an exponential. 

In a finite system, finding the phase transition from the 
crystal state to the SF is not simple, because the finite 
number of spins prevents the system from assuming arbi- 
trary polarizations. This results in a division of the phase 
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diagram into different stripes with fixed integer number 
of up spins, which makes the crystal lobes difficult to 
discern. In infinite systems, however, one expects a step- 
like behavior of the polarization only in crystal phases, 
while in the SF the polarization should change smoothly 
with /i between —1/2 and 1/2. This observation makes 
it possible to extrapolate the border between the crystal 
phases and the SF by the following finite size scaling: At 
fixed and for a given polarization, compute for several 
chain lengths TV the upper and lower limits of the polar- 
ization plateau, /i a (N) and fi^ (N). Then, from a finite 
size scaling of the results, one can extract /i a ,& (N = oo). 
In the SF, where the polarization should be a contin- 
uum, \i a (oo) and fib (oo) should be equal. However, in 
the crystal lobe there will be a finite distance between 
li a (oo) and //& (oo): the width of the lobe for that value 
of 6. A similar procedure is known to work well for the 
estimation of Mott-lobes from finite systems in the Bose- 
Hubbard model |39) . The main difference is that in our 
model there is in principle a lobe for any rational filling 
factor, instead of only for integer filling factors. Panel E 
of Fig. [3] shows the result of this approach for the 1/3- 
filling lobe, using chain lengths of up to 102 spins. The 
cusp structure is typical for one dimensional systems. 



VI. INFINITE TIME EVOLVING BLOCK 
DECIMATION (ITEBD) 

To complement the results of the previous section, we 
study the phase diagram with the iTEBD algorithm [40J , 
in which infinite chains can be directly addressed without 
the need of finite size extrapolations. Treating long-range 
interaction complicates the original formulation of the 
iTEBD algorithm and turns out to lead to convergence 
problems. Instead, we implement a variant of the original 
algorithm in which the translationally invariant character 
of the problem is kept in the state and the interactions, 
thanks to the use of Matrix Product Operators [H] . In 
this way, we can include interaction terms ranging longer 
than nearest neighbors in a much simpler way. 

By comparing the phase diagrams for different range 
of XY interactions, we may study the effect of the long- 
range hopping in the thermodynamic limit . The method 
still requires a truncation of the interacion range to some 
finite order, so that it will not be possible to reproduce 
the power-law decay of correlations within the crystal 
lobes. 

We truncate the ZZ interactions above next-to-nearest- 
neighbor (NNN) interactions and compare the cases of 
NN and up to NNN XX and YY terms. We observe 
that a small bond dimension, D — 10, provides already 
a good approximation to the overall phase diagram. To 
analyze the isolating lobe at 1/3 filling we increase the 
bond dimension to D = 20. 

Figure [5] shows the results for the magnetization per 
particle in the first case, in which XX and YY interactions 
range only to the NN. As in [ST] the phase diagram is 
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Figure 5: Mean magnetization for an infinite chain with NN 
XY interactions and NNN ZZ term. The right panel shows a 
zoom on the isolating lobe at 1/3 filling. 




Figure 6: Mean magnetization for an infinite chain when the 
range of all ZZ and XY interactions is extended to the NNN. 
The right panel shows a zoom on the isolating lobe at 1/3 
filling. 

symmetric. This is also seen in the close up on the filling 
1/3 crystal lobe (right panel of Fig. [5]). 

Figure [6] shows the equivalent phase diagram when 
NNN terms are included in the XY plane. Including 
one more term in the XY interactions already causes a 
clear asymmetry of the phase diagram with respect to the 
change 6* — > — This is clearly visible in the superfluid 
phases (left panel of the figure) , and the zoom on the iso- 
lating lobe at polarization —1/6 shows some asymmetry 
as well. 

The size of the lobe is larger than the result from 
DMRG after finite size extrapolation. We are, however, 
only including the first term further than the NN, and we 
expect that longer range terms correct the exact shape 
of the lobe, as they will give rise to lobes corresponding 
to other filling factors. 

VII. EXACT DIAGONALIZATION 

In this section, we supplement the previous results by 
exact diagonalization (ED) of chains of 18 spins with pe- 
riodic boundary conditions, using the Lanczos method. 
Although ED only allows investigation of very small sys- 
tems, from these one can often infer surprisingly much 
about the behavior of the model at larger sizes. More- 



8 




-0.1 -0.075 -0.05 -0.025 



Figure 7: Fidelity susceptibility \f/W for 18 spins (ED), and 
second derivative of total in-plane magnetization d 2 M x /d9 2 
for 18 spins (ED) and for 60 spins (DMRG). 



over, an experimental implementation with a chain of 
trapped ions will have to start with short chain lengths. 
For such a case ED provides a valuable validation of the 
trapped-ion quantum simulator. 

Figure [7] shows the fidelity susceptibility, 



XF = 



i-iW)iv>(0+A0)>r 



A9 2 



(13) 



for 1 /3 filling. At a second order quantum phase transi- 
tion, xf diverges in the thermodynamic limit. The peaks 
in a finite system are precursors thereof. However, in the 
present system the divergence turns out to be relatively 
weak, which makes it difficult to extrapolate quantitative 
values for the critical point of the transition. The quali- 
tative value, however, is consistent with the DMRG and 
iTEBD analyses. The second derivative of the energy 
has been proposed [H] as a substitute of the fidelity sus- 
ceptibility to detect quantum phase transitions. In the 
present system, however, the energy is very smooth and 
does not give a good indication of the phase transition. 

Instead, we compare in Fig. [7J \f to the sec- 
ond derivative of the total in-plane magneti- 
zation, d 2 M x /d8 2 , where we defined M x = 
Ei*y I (SfSf + SfS*) | / (N(N-l)/2). Here, _ the 
prime indicates that the sum excludes self-correlations. 
We take the absolute value of the correlations in order to 
take into account the different staggered configurations 
occurring for different filling factors. We see that 
d 2 M x /d0 2 peaks at slightly larger absolute values of 
than XFi but the rough locations are consistent with the 
peaks of \f pH] - 

A comparison to the results from DMRG shows a very 
similar qualitative behavior. This means that already the 
very small chains considered here show strong precursors 
of the physics relevant for large chain lengths. 



VIII. CONCLUSION 

To summarize, we have studied a one dimensional sys- 
tem of hard-core bosons where particles interact and can 
tunnel at long distances with an algebraically decaying 
strength. We have discussed how the system can be 
mapped onto a spin Hamiltonian that can be simulated 
experimentally using trapped ions. Unlike other atomic 
quantum simulators (such as dipolar ultracold atoms), 
trapped ions appear to be more flexible in the manip- 
ulation of some parameters — e.g. it is easy to change 
the interactions and hoppings from negative to positive 
simply by detuning the lasers. 

The system we have studied has a rich phase diagram, 
with many insulating phases filled with (possibly long- 
lived) meta-stable states and quasi-long range order in all 
correlations — two prominent features that are induced 
solely by the long-range nature of the interactions. 

The coexistence of diagonal LRO and off-diagonal 
quasi-LRO in the insulating lobes can be seen as a quasi- 
supersolid, the closest analogue to a supersolid that can 
be found in ID. It would be very interesting to study how 
these phases behave in 2D, where the possibility of hav- 
ing off-diagonal LRO could allow the quasi-supersolids to 
become true supersolids. 

Since calculations prove to be quite complicated for 
long-range interactions, one can consider the present 
model as a testbed for a comparison of different numer- 
ical tools. We find that the precision of the standard 
perturbative mean- field theory is insufficient, because it 
does not capture excitations more complicated than sin- 
gle particle excitations. Exact diagonalization, DMRG, 
and iTEBD, however, all deliver a consistent picture. 
Among these, we find that DMRG is the method of choice 
for the model studied. However, it is restricted to finite 
systems. While iTEBD is less accurate and restricted to 
shorter range of interactions, it allows to investigate the 
thermodynamic limit directly. Exact diagonalization, fi- 
nally, while being relevant only for very small systems, 
can still serve for benchmarking other less accurate com- 
putational methods and near-future experiments. 

In the future, it would be interesting to pursue investi- 
gations of other properties of the system that are affected 
by the long-range nature of the interactions, such as its 
response to excitations or the dynamics of its correla- 
tions. 
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